Effective approach to the Nagaoka regime of the two dimensional t—J model 
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We argue that the t—J model and the recently proposed Ising version of this model give the same 
physical picture of the Nagaoka regime for J ft <C 1. In particular, both models are shown to give 
compatible results for a single Nagaoka polaron as well as for a Nagaoka bipolaron. When compared 
to the standard t—J or t-J z models, the Ising version allows for a numerical analysis on much larger 
clusters by means of classical Monte Carlo simulations. Taking the advantage of this fact, we study 
the low doping regime of t—J model for J/t <C 1 and show that the ground state exhibits phase 
separation into hole-rich ferromagnetic (FM) and hole-depleted antiferromagnetic (AFM) regions. 
This picture holds true up to a threshold concentration of holes, 8 < 5 t ~ 0.44-y/ J/t. Analytical 
calculations show that S t = Jj/2-irt. 



I. INTRODUCTION 

It is by and large recognized that the key proper- 
ties of high-temperature superconducting materials can 
be explained with the help of the two-dimensional t—J 
modelPEI Therefore, the major interest is focused on 
this model in the regime that seems to be relevant to 
high-temperature superconductivity. Since the Hubbard 
on-site repulsion U in cuprates is usually assumed to be 
an order of magnitude bigger than the hopping integral 
t, the relevant AFM coupling between nearest neighbor 
sites J = At 2 /U is about tenths of t. 

Much less is known about the t-J model in the small- 
J limit, i.e., the large-// limit of the Hubbard model. 
At half filling the infinite-?/ Hubbard model has AFM 
ground state. However, in 1965 Nagaoka proved a 
theorem- which states that when exactly one hole is in- 
troduced the ground state becomes FM. In the infinite-// 
limit the ground state of the half filled Hubbard model 
is macroscopically degenerate. When a single hole is in- 
troduced this degeneracy is lifted, since it is energetically 
favorable for the hole to move in a background of fully 
aligned spins. A simple proof of Nagaoka's theorem was 
later given by Tasaki,^ who also showed that additional 
density-dependent interactions do not alter this result. 

The Nagaoka theorem is one of very few rigorous re- 
sults concerning strongly correlated electronic systems. 
However, it does not say anything about the case of a 
finite density of holes as well as finite AFM interaction. 
The question of a character of the leading instability of 



the Nagaoka state with respect to the AFM exchange 
term or finite density of holes has attracted much in- 
terest. For finite J (J > 0) and/or finite doping the 
ground state is determined by the competition between 
antiferromagnetism favored by the exchange interaction 
and Nagaoka's ferromagnetism favored by the kinetic en- 
ergy. This competition presumably drives the system 
into two phases, a hole-rich FM region where the kinetic 
energy is minimized and a region with localized electrons 
characterized by AFM order. Generally, the t-J model 
may display different types of phase separation depend- 
ing on the value of Jit and the actual state of affairs is 
far from being clearP^It was demonstrated in Refs. [B]- 
IHI that phase separation takes place in the t—J model for 
all values of J. Other authors^"^ g nc j phase separation 
only for large J. It is very difficult to establish unambigu- 
ously the presence or absence of phase separation in the 
small-J limit of the t-J model. The main reason is that 
a high-accuracy, unbiased calculation of the ground state 
energy as a function of the hole density is required to as- 
sess the competition between the interaction and kinetic 
energies. The spatial inhomogeneity in the phase sepa- 
rated state makes analytical approaches rather involved. 
On the other hand, since the FM bubbles are relatively 
largej^ it is difficult to apply numerical approaches like 
the quantum Monte Carlo method or exact diagonaliza- 
tion. In this situation, a computational method that is 
able to tackle a system sufficiently large to describe the 
spatially separated state in the small-J limit of the t-J 
model is required. 
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In this paper, we demonstrate that Monte Carlo simu- 
lations for the recently proposed Ising version of the t—J 
modeP^ provide sound and reliable results in this limiting 
case. We employ this approach to study the formation 
of a bubble of the FM phase when holes are introduced 
into an AFM background. Namely, we investigate the so- 
called Nagaoka polaron which sets in for vanishing doping 
and J/t <C lP^ The Nagaoka polaron represents an in- 
termediate state between the homogeneous FM Nagaoka 
ground state for J — and the standard spin polaron 
for J > 0.2. Numerical studies of the Nagaoka polaron 
are difficult because of its large spatial dimensions: for 
J — > its radius diverges as J~ 1//4 . 

In this paper we formulate an effective description of 
the Nagaoka regime, which is based on the recently pro- 
posed Ising version of the t-J modelP-^In the subsequent 
sections we recall the basic properties of this model and 
explain why it gives the same physical picture of the Na- 
gaoka regime as the standard isotropic t-J model. These 
qualitative arguments are accompanied by quantitative 
comparison with the numerical results for the isotropic 
model. For the reader convenience we first recall density 
matrix renormalization group (DMRG) studied^ on the 
Nagaoka polaron. We then present new results for the 
Nagaoka bipolaron, studied in the t-J model by means 
of exact diagonalization in the limited functional space 
(EDLFS)P 

Numerical calculations within the Ising version are by 
far less demanding hence much bigger clusters and/or 
much larger doping become accessible. After having suc- 
cessfully tested the single and two-holes cases, we in- 
vestigate the ground state properties of the t-J model 
(J/t <C 1) doped with several holes. In particular, we 
show that all holes are confined in a single FM polaron. 
We discuss how its energy and spatial dimensions depend 
on the number of holes. For low hole densities, our re- 
sults provide an evidence that the leading instability of 
the FM Nagaoka state is a phase separation rather that 
a single spin flip. 

II. THE ISING t-J MODEL 

We start with the t-J Hamiltonian on a square lattice^ 

H t -j = -J2 + J zZ [QiQi - I^j) > (!) 

ijtr (ij) 

where c ia = Pci a P = 0^(1 — n^-a) is the projected elec- 
tron operator (to exclude the on-site double occupancy) , 

Qi = J2cr,a> c \a T <y<j'C iiy >, t 2 = 3/4, is the electron spin 
operator and n, = PriiP = + — 2n^n^ is the 
projected electron number operator. Hamiltonian ([I]) 
contains a kinetic term with the hopping integrals tij 
and a potential J describing the strength of the nearest 
neighbor spin exchange interaction. At every lattice site 
the Gutzwiller projection operator P = JX^l — n it7 n i _ a .) 
projects out the doubly occupied states. Physically this 



modification of the original Hilbert space results in strong 
electron correlation effects. 

At the low-energy scale of order J (<C t) , it is reason- 
able to consider the background spin configuration to be 
nearly frozen with respect to the hole dynamics. In this 
case the properties of the low-energy quasiparticle excita- 
tions in the t-J model are at least qualitatively similar to 
those in the anisotropic t-J z model, in which the spin-flip 
part of the Heisenberg interaction is dropped: 

H t-J- = ~zZ tij^rCir + J, \QiQ Z j - > ( 2 ) 

ija (ij) V J 

The global continuous spin SU(2) symmetry of the t-J 
model now reduces to the global discrete Z2 symmetry 
of the t-J z model. Although QfQj interaction possesses 
discrete Z2 symmetry, the original SU(2) symmetry of all 
other terms of the Hamiltonian is preserved. Therefore, 
the symmetry of the t-J z model depends on whether J z 
is zero or finite. Namely, for J z = the SU(2) symmetry 
is restored again. Although this model is more amenable 
to numerical calculations, again only rather small lattice 
clusters are allowed. 

One may hope that the full Ising version of the t- 
J model in which the i-term also possesses the global 
discrete Z2 spin symmetry results in a more tractable 
though still nontrivial model. It by definition has the 
global discrete Zi symmetry, regardless of the values of 
the incoming parameters. This implies that the result- 
ing system can be thought of as a doped classical Ising 
model. However, it is not clear how such a model can 
be derived directly from ([!]), since the projected electron 
operators Ci„ transform themselves in the fundamental 
representation of SU(2). To overcome this problem, the 
recently proposed spin-dopon representation of the pro- 
jected electron operators can be usedP^ 

The idea behind that approach is to assign fermion 
operators di a to doped carriers (holes, for example) 
rather than to the lattice electrons. The enlarged on- 
site Hilbert space is spanned by the state vectors \cra), 
with a =ft,J| labeling the 2D Hilbert space of the lat- 
tice spin, Si, and a = O,t,4ot4- labeling the AD on-site 
dopon Hilbert space. The physical space consists of the 
spin-up I 0)j spin-down | -IJ- 0)j and spinless vacancy 
(I lU)i - I Wi)/V2 statesPl The remaining unphysical 
states are removed by the constrainPS 

S t M t + ^nf = 0, (3) 

where Mi — a' ( $i a T <ja'd'io l stands for the dopon spin 
operator so that 

Qi = Si + Mi. (4) 

The physical electron operator Ci a is then expressed in 
terms of the spin and dopon operators projected onto 
the physical subspace singled out by Eq. 
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In view of the relation (S a ) 2 — 1/4, the constraint ^ 
can equivalently be written in the form 



S?M? 



(S?) 



0. 



(5) 



Within the full Ising t-J model, one should have Q~ = 
Qf ± Q\ = 0. In view of Eq. (4|), this requires Sf = 



Mf 



0. To explicitly derive t 



re Ising t-J model, we 



therefore project the dopon operators onto the Hilbcrt 
space singled out by the local constraint 



SfMf 



1 



0. 



(6) 



which can be thought of as an "Ising" form of Eq.Q. It 
represents the Z2 singlet under the Q\ — > ±Qf transfor- 
mations e Z2 C SU(2). The physical electron projected 
operators reduce to the Z2 spinors: 



2 s * 



if 



(7) 



'vi = n ";-|A-; = I - + t v- j (8) 

with the projection operator Tf h = 1 - (2Sf Mf + nf /2). 
It can readily be checked that 



1 



QJ = (QT) 



0. 



as desired: the transverse components of the electron spin 
operators no longer appear in the theory. 

The underlying onsite Hilbert space rearranges itself in 
the following way. The operators §7§ act on the Hilbert 
space Hi = {| JJ-,0), | -D-,t)}- These operators do not mix 
up any other states. Operator Cjj, destroys the spin- 
down electron and creates a vacancy. This vacancy is 
described by the state | -IJ-,t)- The similar considera- 
tion holds for the operators. Now, however, the va- 
cancy is described by the state | f|",|). Those two va- 
cancy states are related by the Z2 transformation. The 
operator (Qf) 2 = \(l — nf) produces zero upon acting 
on the both. The physical Hilbert state is therefore a 



Hi 



H^. Under the Z 2 transforma- 



direct sum H p h 
tion Sf -t —Sf) we get H\ O Hi, which re- 

sults in H p h —> H p h- In the isotropic t-J model these 
two 2D spaces merge into a 3D SU(2) invariant phys- 
ical space, where the vacancy is just an antisymmet- 
ric linear combination given by the SU(2) spin singlet, 
(I 1tl)i ~ I -JJ-T)-i) / "s/2- The symmetric combination splits 
off, since it represents an unphysical spin-triplet state. 

As a result, one arrives at the representation ^ in 
which, however, the electron projection operators are 
given by Eqs. (7 [8|. All the parts of this Hamilto- 
nian possess the global discrete Z2 symmetry whereas 



the global continuous SU(2) symmetry is completely lost. 
Close to half-filling, the Ising t-J Hamiltonian takes on 
the formpS 



H, 



Ising 
t-J 



ija 



<ij> 



5Z QZ 
i " 1" 



+S?M* 



S?Mf 



(9) 



which is to be accompanied by the constraint ([6]). The 
magnetic M?M? term has been dropped as being small 
of order 5 2 in the limit S := (rid) *C 1. 

In practical calculations, we find convenient to imple- 
ment the constraint ^ with the help of a Lagrange mul- 
tiplier. Since SfMf + nf/4 > 0, the global Lagrange 
multiplier 



SfMf 



+ 7 n 
4 % 



(10) 



enforces the constraint ^ locally. The unphysical occu- 
pancy of an arbitrary site would enhance the total energy 
by A —¥ +00. Therefore, all unphysical on-site states are 
automatically eliminated, so that the local constraint is 
taken into account rigorously. 

The main difference between the t—J z and Ising t-J 
models originates from the different symmetries of the 
hopping terms as discussed inP^ The one-hole energy 
obtained for the isotropic t-J model is shown to be in 
between the results obtained for the t-J z and Ising mod- 
els. In the regime J <C t, the differences between t-J 
and t-J z models are comparable to those between t-J 
and Ising models. 



A. Monte Carlo approach 



The Hamiltonian (JoJ) together with the constraint ( 10 ) 
describe a system that contains both classical (Sf) and 
quantum (di a ) degrees of freedom. However, there is no 
direct interaction between the quantum particles. There- 
fore, the Ising t-J model is closely related to a (multi- 
component) Falicov-Kimball model and we can apply an 
efficient Monte Carlo (MC) approach derived exclusively 
for the latter model p2 This numerical approach can be 
applied neither to the standard t-J model nor to the t-J z 
one. In the latter case one should use the SU(2) invariant 
constraint ^ which involves S ± operators. 

The classical MC method has already been applied to 
the Ising t-J model in order to study dynamics of holes 
and destruction of the AFM order with increasing hole 
concentration.^ Here, however, we are not interested 
in the thermodynamics of the system, but rather in its 
ground state properties. Therefore, the Metropolis algo- 
rithm is used for simulated annealing.^ And once again, 
since simulations are performed by random walk through 
the configuration space of the classical variables, there is 
no need for quantum annealing and relatively large lat- 
tices can be studied. The size of the lattice for a given 
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J /t and a given number of holes is adjusted so that the 
size of the polaron be always significantly smaller. Since 
the holes can propagate only within the FM region, the 
finite size effects become negligible. Most of the calcula- 
tions are carried out on 20x20 and 30x30 lattices, but in 
some cases larger clusters, up to 50x50 are necessary as 
well. From now on, a value of the nearest neighbor hop- 



ping amplitude is used as the energy unit (U 



t = 1 for 



i and j being nearest neighbors and f,-j = otherwise) 



III. TEST OF A SINGLE POLARON PROBLEM 

Neglecting the spin-flip term is generally a crude ap- 
proximation, hence the standard and the Ising t—J mod- 
els describe very different systems. However, we will show 
that for vanishing doping and small J both approaches 
give rise to the same physical picture of the Nagaoka po- 
laron. 

The density matrix renormalization group (DMRG) 
studies®! of the 2D t—J model have shown that for 
J < 0.03 the Nagaoka polaron is indeed stable. Its size 
and energy can be determined by balancing the kinetic 
energy of a hole freely propagating within a FM bubble 
against the magnetic energy of the FM bubble relative to 
the energy of the Neel state. Minimizing the sum of these 
two energies one easily finds expressions for the radius R 
and energy E of the Nagaoka polaron (see Ref. [T8]) : 



R~ 1.12J 



-1/4 



E ~ -4 + 9.2VJ, 



(11) 



which for J < 0.03 accurately fit the DMRG data. 

The Ising t-J model displays essentially the same 
physics. The constraint ( 10 ) allows for propagation of 



holes only within static FM bubbles where S? = — M? 
and formation of these bubbles takes place at the expense 
of the exchange interaction, which favors AFM alignment 
of Sf . In both models the magnetic energy of a bubble is 
qualitatively similar and is proportional to J multiplied 
by the number of FM bonds. Quantitative differences 
should arise from different energies per bond of the Neel 
ground state of the undoped systems. Figs. [l]and[2]show 
quantitative comparison of both models. Here, we show 
the ground state properties of a single Nagaoka polaron 
in the Ising t—J model for a 50 x 50 system with peri- 
odic boundary conditions and A = 300. The radius and 
energy of the Nagaoka polaron have been compared with 
expressions (111 as well as with the bare DMRG data for 
the SU(2) t-J model taken from Ref. HI The overall 
agreement is clearly visible. 

To conclude this section, we notice that the Ising t-J 
model provides a simple and reasonably accurate descrip- 
tion of a single Nagaoka polarons for J <C 1 simply be- 
cause the spin-flip term becomes irrelevant inside a suf- 
ficiently large FM bubble which confines the movement 
of holes. 
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FIG. 1: (Color online) J-dependence of the radius of the Na- 
gaoka polaron formed by a single hole. Points show results 
from the MC calculations for the Ising t-J model, continuous 
line is the power-law fit a nd t he dashed line shows the de- 
pendence described by Eq. (111. Inset shows the same results 
but on the log-log scale. 
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FIG. 2: (Color online) The same as in Fig. [T] but for the 
energy of a single Nagaoka polaron relatively to the energy 
of the homogeneous Neel state. The energy is compared with 
DMRG results for the t-J model taken from Fig. 4 of Ref. 
1181 The other lines correspond to Eq. (Ill for the isotropic 



t—J (dot-dashed violet line) and a similar expression for the 
t—J z model from Ref. 18 (dashed green line). The continuous 
line shows a power-law fit to the present results. 



IV. THE NAGAOKA BIPOLARON PROBLEM 

In the preceding section we have successfully tested the 
case of a single carrier doped into Mott insulator. Since 
our aim is to study a system at finite doping such a test 
might still be insufficient. Therefore, we investigate also 
the system doped with two holes. On the one hand it 
is a first nontrivial step towards understanding the spa- 
tial hole distribution in doped systems (uniform versus 
inhomogeneous). On the other hand, due to large spatial 
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dimensions of the Nagaoka polaron, it is already a chal- 
lenging problem for fully quantum numerical approaches. 

To investigate the Nagaoka bipolaron in the isotropic 
t-J model, we employ the EDLFS method. This method 
describes properties of carrier/carriers doped into a pla- 
nar ordered antiferromagnetP^ One starts from a trans- 
lationally invariant state of two carriers in the Neel back- 
ground 



6 o) P = E(- 1 ) M(T)c « c -vl Nee1 )' 



(12) 



where the sum runs over two nearest neighbors to site 
and M(-y) sets the appropriate sign to generate p x (y)- 
wave symmetry. The kinetic part Hk as well as the off- 
diagonal spin-flip part Hj of the Hamiltonian are 
applied up to Nh times generating the basis vectors: 

{\^)} = [H k + Hj] n »\<j> ) p , n h = 0,...,N h . (13) 

Then, the ground state l^o) is calculated within the lim- 
ited functional space by means of the Lanczos method. 
The advantage of EDLFS over the standard exact diag- 
onalization (ED) approach follows from systematic gen- 
eration of selected states which contain spin excitations 
in the vicinity of the carriers. It enables investigation of 
much larger systems, which is particularly important in 
the Nagaoka regime. We apply this method and study 
the average distance D between two holes in the t-J 
model. 

From the construction of EDLFS follows that Nh de- 
termines maximum distance between holes and spin exci- 
tations as well as the maximum accessible D. Therefore, 
all characteristic length-scales of investigated problem 
should be smaller than a certain £(Nh). Since the suc- 
cessive application of the nearest neighbor hopping [see 
Eq. ( 13 )] closely resembles the random walking process 



we expect that £,(Nh) oc \/Nh. The numerical complex- 
ity of the two-holes problem allows us to study up 
to 12. Therefore, for small J we carry out a finite size 
scaling with respect to the size of the Hilbert space and 
study D(J) = D(J,Nh — > oo). The finite size effects 
should vanish when £(Nh)/D(J) > 1 or equivalently 
7V/ t /D 2 (J) ^> 1. We have found that this vanishing can 
universally be described by the Gaussian: 



D(J)-D(J,N h ) = aexpi-l3 



N h 



D 2 (J) 



(14) 



where parameters a and (3 are independent of J. Note 
that for a given J fitting of D( J, Nh) involves only a single 
free parameter D{ J) which simultaneously represents the 
extrapolated distance between two holes. Fig. [3] shows 
the extrapolation, while the resulting D{J) is shown in 
Fig. [4] Finally, the quality and the universality of the 
finite-size scaling is directly shown in Fig. [5] One can 
clearly see that all data for J = 0.04, 0.14 merge 
into a single curve, which strongly supports the assump- 
tions behind Eq. ( 14 ) . This also demonstrates a self- 
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FIG. 3: (Color online) Distance between two holes D in the 
isotropic t-J obtained in a functional Hilbert space gener- 
ated for finite Nh- Numerical results (points) are fitted and 
extrapolated according to Eq. (14 1. 
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FIG. 4: (Color online) Extrapolated distance between two 
holes D(J) for the isotropic t-J model. 



These results provide solution of the long-standing 
open problem, i.e., whether two hole s in t he t-J model 
form a bound state in the small J limit JM2H Since D( J) > 
3.5 for J < 0.10, the problem can hardly be solved by ex- 
act diagonalization on 32-site cluster^ where the max- 
imal possible distance between two holes is four lattice 
spacings. Note that the symmetry of the bound state 
is yj-wave for J < 0.15. Therefore, the results using the 
EDLFS method indicate that two holes in the single band 
t—J model are always bound, which may not necessarily 
be the case in more general models describing the C11O2 
planeP 
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FIG. 5: (Color online) The same data as in Fig. [3] but 
rescaled according to Eq. 1 14 1. 



The most important result shown in Fig. [4] concerns 
the power-law dependence: 



D(J) ~ 1.93J 



-0.27 



(15) 



hence D(J) is roughly proportional to the radius of a 
single-hole polaron R{J) [see Eq. (11)]. Already this 



result suggests that the linear dimensions of polaron and 
bipolaron are determined by the same mechanisms, which 
are properly captured by the Ising t-J model. In order 
to show that this expectation holds true we have calcu- 
lated D in the Ising version of the t-J model. Results 
are shown as (red) stars in Fig [6} In the same figure re- 
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FIG. 6: (Color online) Comparison of the distance between 
two holes in the t-J (blue circles) and Ising t—J (red stars) 
models. "Ising t—J rescaled" means that all values are multi- 
plied by 1.33. Note that the rescaled values can be described 
by Eq. ( 15 1, see also the inset where compounded D( J) data 
are presented on the log - log scale. 



suits for the t-J model are presented as (blue) dots. One 
can clearly see that the distance between holes D in the 
Ising t-J model increases with decreasing J slower than 
in the full t-J model. This result is not a surprise: in the 
isotropic t-J model (as well as in the t—J z one) a hole 
is able to enter the AFM surrounding of the FM bubble. 
Since the boundary of the FM bubble is unpenetrable 
in the Ising t—J model, we expect a smaller average dis- 
tance between holes in the latter case. However, it turns 
out that the difference has only a quantitative charac- 
ter, i.e., only the coefficient in Eq. |l5]) is approx. 33% 
smaller. This agreement shows that these two methods, 
i.e., the EDLFS method for the full t-J model and the 
MC method for the Ising t-J model are complementary 
in a sense: The applicability of the EDLFS method is lim- 
ited by the maximum size of the Hilbert space, and since 
the distance between holes increases with decreasing J, 
this method cannot be used when J is too small. On the 
other hand, the importance of the spin-flip term in the t— 
J Hamiltonian diminishes with decreasing J and the ap- 
proximation that leads to the Ising t—J model becomes 
more reliable in the region where the EDLFS method 
cannot be applied any longer. The main advantage of 
the Ising t—J model is that it can be studied within the 
framework of the classical MC method on clusters suf- 
ficiently large to describe large polarons that emerge at 
small J. 



V. FINITE HOLE DENSITY 

Up to this point we have been analyzing one and 
two holes in the whole system. Since the size of the 
(bi)polaron and its energy do not depend on the size 
of the lattice [provided the lattice is significantly larger 
than the (bi)polaron size], these results effectively de- 
scribe the case of the vanishing density of holes. Then, 
the important question arises as to how the ground state 
of the Ising t-J model evolves when the number of holes 
increases. Possible scenarios include phase separation or 
homogeneous distribution of holes. It is also possible that 
a single polaron becomes unstable at some critical value 
of the hole number giving way to smaller polarons. 

In order to study this problem we calculate the total 
energy of the system as a function of the number of holes 
E(N). Convex E(N) for some N indicates that it is 
energetically favorable to split a FM bubble with N holes 
into smaller bubbles with M < N holes, provided that 
E(M) is concave. If E(N) is convex for arbitrary N, 
holes will not form polarons with more than one hole. 
On the other hand, if E(N) is concave for arbitrary N, 
all holes introduced into the system will gather in a single 
FM region. In other words, the phase separation into a 
hole-rich FM region and an AFM region without holes 
takes place. 

Accurate determination of E(N) is generally a diffi- 
cult task. For N = 1 the shape of the polaron is almost 
circular, but for higher N the geometry becomes nontriv- 
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ial. Fig. [7] shows the Nagaoka polarons and correspond- 
ing hole wave functions for N = 2, . . . , 5. The shape 
of the polaron follows directly from the spatial structure 
of the occupied orbitals. The diagonal orientation of al- 
most rectangular polarons minimizes the magnetic en- 
ergy along the line between FM and AFM regions. 

With increasing N the size of the FM bubble increases, 
so a large lattice is necessary to avoid the finite size ef- 
fects. The advantage of the Ising version of the t-J model 
becomes evident, since it can be reliably investigated on 
lattices much larger than those accessible to the fully 
quantum methods like quantum MC, exact diagonaliza- 
tion or even EDLFS. Using larger lattices we study po- 
larons containing up to 10 holes. In Fig. [8] we show 
the polaron energy E(N) as a function of the number 
of holes N for J = 0.01. Studying other values of J 
from 0.01 to 0.1 (not shown) we fitted the energy with 
a function E(N) = aN + by/N. In Sec. V.A. we justify 
such a form of E(N). We have found that in all cases b 
is positive, which means that E(N) is concave. This in 
turn implies that the Ising version of t-J model displays 
a phase separation for all those values of J at which it is 
still equivalent to the isotropic t-J model. 

An important question concerns the fraction of the sys- 
tem occupied by each of the magnetic phases. It can be 
answered by comparing the size of the polaron to the size 
of the whole system. Fig. [9] shows the relation between 
the fraction of the lattice sites with ferromagnetically 
aligned spins N p /N s and the density of holes S = N/N s . 
N p is the number of lattice sites in the polaron and N s is 
the size of the lattice. This dependence can be fitted by 
a linear function, which extrapolated to N p /N s = 1 gives 
the threshold value of the hole density St ■ If the concen- 
tration of holes is close but still smaller than S t , most of 
the system is occupied by the FM phase, while the rest 
forms an AFM island (or islands) . Finally, for concentra- 
tions larger than St the whole system would be in a fully 
polarized state. However, the latter regime is probably 
not accessible by the present approach. In the Nagaoka 
regime, the Ising and isotropic t-J models give the same 
results because the physics of the Nagaoka regime is de- 
termined by the competition between the magnetic and 
kinetic energies. However, as soon as the FM bubble cov- 
ers the whole system other mechanisms come into play, 
e.g., a direct hole-hole interaction and/or interference of 
the carriers paths around loops. Fig. [10] shows St as 
a function of J (the point St — 0, J — is a result of 
the Nagaoka theorem). The obtained square-root depen- 
dence between both quantities follows immediately from 
Eq. ( 11 ) and the proportionality N p oc N shown in Fig. 
1 



A. Analytical approach for finite doping 

We consider the FM polaron of size N p with N doped 
holes that can be treated as spinless noninteracting 
fermions. The FM polaron is furthermore placed in the 



hole-depleted Neel spin background. We furthermore 
consider the limit of small hole-density that allows for 
quadratic expansion of the single particle kinetic energy: 



^kin( fc ) = -2( c °s k x + cos ky) 



-4 + fc 2 . 



(16) 



We obtain the kinetic energy of N holes by integrating 
Eq. Q up to k F = 2^/itN/N p 



E kin = -4N + 2itN 2 /N p . (17) 
We proceed by writing the total energy E(N) as: 



E(N) = E kh 



E. 



spin j 



E(N) 



-4N + 2nN 2 /N p + N p J, 



(18) 



where the last term represents the magnetic energy of 
the FM polaron relative to the energy of the Neel state. 
After the minimization 8E/dN p = we obtain 



Np 



S, 



(19) 



representing the ratio between the polaron size N p and 
the total size of the system N s as a function of hole doping 
S. The comparison of Eq. ( 19 ) with numerical data is 
shown in Fig. [9] From Eq. ( 19 1 we obtain as well the 
critical doping for the transition to the FM state 



S t = 



(20) 



shown along the numerical results in Fig. [TO] Notice, that 
Eq. ( 20 ) agrees with that derived within the semiclassi- 
cal calculations of the 2D isotropic t-J modeP^ which 
suggest that at small hole concentration and rather weak 
AFM coupling the FM Nagaoka state becomes unstable 
towards a creation of an AFM bubble. This agreement is 
quite natural, since the spins are considered to be frozen 
in both the classical large-spin limit of the isotropic t-J 
model and its full Ising version. Finally, we obtain the 
total energy of the system 



E(N) 



-At 1 




(21) 



The comparison with numerical results is shown in 
Fig. [8j Here we have neglected the effects along the 
line, separating the FM polaron from the Neel spin back- 
ground as well as the dependence of the kinetic energy 
on the shape of the bubble. As a results only the linear 
term in E(N) is reproduced. Since the phase separa- 
tion is determined by the nonlinear part of E(N), these 
effects give rise to small, nevertheless important correc- 
tions. The former one is proportional to the length of the 
borderline (oc y/~N) and it was the reason for the choice 
of the fitting function in Fig. [8] 
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FIG. 7: (Color online) Leftmost column: shapes of the Nagaoka polarons including from 2 to 5 holes. Filled circles indicate 
spin-up lattice sites and empty circles spin-down lattice sites. The rest of the panels show wave functions of the holes. In all 
cases J=0.03 was assumed. 
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FIG. 8: (Color online) Energy E(N) of a polaron containing 
N holes relatively to the energy of the homogeneous Neel 
state (points) for J = O.Olt. The full line represents a fit to 
the numerical data as given in the legend, the dashed line 
represents analytical result in Eq. (211. 
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FIG. 9: (Color online) Fraction of the lattice occupied by the 
FM polaron as a function of the concentration of holes for 
different values of J. The straight lines represent analytical 



results given in Eq. ( 19 1 



VI. SUMMARY 

The main difficulty in analyzing the t-J model in the 
small- J limit is that a large size of the lattice is required 
to correctly describe the dynamics of holes. This require- 
ment significantly restricts the applicability of numerical 
approaches like the quantum MC or exact diagonaliza- 
tion method. 

In the small- J regime, however, the holes are con- 
fined in a FM polaron, so that the spin-flip processes 



FIG. 10: (Color online) Critical value of the concentration 
of holes above which the whole system is in a fully polarized 
FM state. The point at J = is added as a result of the Na- 
gaoka theorem. The dashed line represents analytical result, 
Eq. pu). 



are strongly reduced. This justifies the applicability of 
the Ising version of the t-J model to study the small-J 
limit of the original t-J model. 

For small but finite values of J, our results for one and 
two holes are in a good agreement with those obtained 
within the fully quantum approaches (DMRG, EDLFS). 
However, we are able to extend our calculations to the 
regimes of smaller J and larger number of holes inacces- 
sible by the former methods. We show that it is ener- 
getically favorable for the system to segregate into the 
FM hole-rich phase and hole-depleted AFM phase. The 
size (surface) of the FM bubble depends linearly on the 
number of holes while its dependence on J is given by the 
square-root function. With increasing concentration of 
holes and/or with decreasing J the size of the FM polaron 
increases and eventually for S t — 0.44yJ it occupies the 
whole lattice. Our numerical results thus suggest that 
Nagaoka state breaks down by forming an AFM bubble. 
This observation fully agrees with a conjecture discussed 
earlier within the isotropic t-J model™ We, however, 
expect that the results obtained for isotropic and Ising 
t-J models start to deviate from each other when dop- 
ing becomes larger that the threshold density St even for 
J<1. 

A rather simple analytic treatment of the holes doped 
in the FM polaron that is furthermore placed in a Neel, 
hole-depleted spin background, leads to a good agree- 
ment with numerical data. Among other results, it 
provides a simple expression for the threshold density 
5 t = \J J/2ir. The theory reproduces only the linear 
dependence of the total energy on the number of holes 
and does not provide information on whether the system 
phase separates. Corrections (e.g., due to line contribu- 
tions) are expected to give rise to a positive -\/lV term in 
the total energy, as obtained from the numerical data. 
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